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Spatio-temporal imaging of light propagation is very important in photonics because 
it provides the most direct tool available to study the interaction between light and its 
host environment. Sub-ps time resolution is needed to investigate the hne and complex 
structural features that characterize disordered and heterogeneous structures, which are 
responsible for a rich array of transport physics that have not yet been fully explored. 
A newly developed wide-held imaging system enables us to present a spatio-temporal 
study on light transport in various disordered media, revealing properties that could not be 
properly assessed using standard techniques. By extending our investigation to an almost 
transparent membrane, a conhguration that has been difficult to characterize until now, 
we unveil the peculiar physics exhibited by such thin scattering systems with transport 
features that go beyond mainstream diffusion modeling, despite the occurrence of multiple 
scattering. 


Introduction 

Light represents a useful and versatile probe to study our world. By traeking its evolution in 
the domains of time, spaee, frequeney and reeiproeal spaee, we are able to extraet a wealth of 
information regarding the surrounding environment. Further insight is often gained when we 
are teehnieally able to traek this evolution in more than one sueh domain simultaneously, as 
demonstrated by a number of reeent results from diverse researeh fields^ Of the different 
approaehes, spatio-temporally resolved teehniques undoubtedly offer the most straightforward 

‘Corresponding author 

^Current affiliation: Laboratoire Kastler Brossel, UMR 8552, CNRS, Ecole Normals Superieure, Universite 
Pierre et Marie Curie, College de France, 24 rue Lhomond, 75005 Paris, France 
^Current position: Gestione SILO Sri, via di Castelpulci 14/d, 50010 Scandicci (FI), Italy 


1 



approach, as underlined by the continuously growing interest in fast photography^ and 
imagingapplications. Further development of similar techniques is key to our ability to in¬ 
vestigate increasingly more complex media, such as disordered or heterogeneous media, which 
are ubiquitous in most helds of study, from atmospheric physics to biology and cultural heritage 
preservation. In this context, one main task of optical investigations is to characterize the light 
transport process determined by multiple scattering. This characterization is typically achieved 
through repetitive measurements to retrieve meaningful microscopic transport parameters, such 
as the transport mean free path k and the absorption length 4, which are directly related to 
the properties of the investigated medium (that is, the internal structure, composition and/or 
density), and are key parameters to consider from both the fundamental and the application 
points of view. 

Here, we experimentally combine wide-held spatial imaging and sub-picosecond (ps) time 
resolution and show how the acquired multi-dimensional dataset allows us to characterize 
previously inaccessible scattering systems and to improve the accuracy of measurements in 
standard situations. The ultrafast time resolution needed for such applications represents a 
major experimental challenge. Established fast imaging techniques involve electronic-based 
devices such as streak cameras or sequential scanning over a given sample in a pump-probe 
conhguration. Despite providing very accurate results, many of these techniques are limited by 
low-time resolutions^^ or a long acquisition time^^. Conversely, wide-held images at sub-ps 
time scales are easily obtained with optical gating techniques; however, their limited accuracy 
has limited their use to qualitative pattern/shadowgram recognition rather than to absolute 
intensity measurements across the entire held of view. In other words, optical gating has been 
previously exploited as a convenient method to see through a turbid medium obliterating the 
complex light transport arising from its structure, not to study it directly. 

In this study, we investigate light transport focusing on three classes of disordered media that 
are extremely relevant for practical applications and yet still very difficult to characterize accu¬ 
rately, namely samples with limited thickness, large-scale heterogeneities and semi-transparent 
membranes. In each of these cases, transport properties are retrieved by exploiting, for the hrst 
time, an all-optical gating setup, which we developed to offer both wide-held acquisition and 
high quantitative accuracy on a sub-ps time scale. Our spatio-temporal investigation allows 
us to grasp the evolution of light transport in its entirety, consequently revealing effects that 
were either overlooked or misinterpreted to date, including a peculiar regime occurring in 
ultra-thin disordered media that becomes apparent only when the spatial and temporal domains 
are studied simultaneously and would not be detectable with other state-of-the-art techniques. 
By comparing the outcome of our analysis with that of the mainstream evaluation model based 
on the diffusive approximation (DA), we describe how single-domain investigations performed 
to date are susceptible to particularly deceptive artifacts that can be clearly appreciated in a 
multi-domain dataset. At the same time, as widely suggested in the literature, we discuss how 
much more robust interpretations, even within the diffusive framework, can be given in terms 
of transverse transport, to which we gain direct access through our technique. 
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Figure 1: Experimental setup and spatio-temporal resolution, a Two linearly polarized fs-pulses 
are provided by a Ti:Sa laser (probe in red) and a synchronously pumped optical parametric oscillator 
(gate in gray). The probe beam interacts with a slab sample, which is imaged on a BBO crystal with 
a 4/ lens configuration. A small aperture (FF) allows the k vectors to be isotropically filtered, while 
a dichroic mirror (M2) guarantees signal and gate spatial superposition. A long-pass (F1) and a 
bandpass filter (F2) block the visible and unconverted light, respectively, b A typical cross-correlation 
time-resolved measurement acquired with a PMT. A fit with the cross-convolution of two sech^ 
pulses returns a FWFIM of 174fs. c The spatial information carried by the pulse is retrieved using a 
CCD camera, whose resolution is influenced by the sum-frequency process, d A geometric sketch 
qualitatively showing how the phase-matching condition results in different angular acceptances in 
the extraordinary and ordinary directions. A slightly misaligned kg and a largely misaligned might 
give the same A6, therefore explaining the sharper resolution in the y direction. 

Materials and Methods 

Experimental setup The optical scheme employed to resolve the time and space propagation 
of light in complex media is based on the cross-correlation gating techniqueshown 
schematically shown in Fig. la. Two synchronous, collinear probe and gate pulses at different 
wavelengths impinge on a 2 mm thick p-Barium borate (BBO) crystal. When the two beams 
overlap both spatially and temporally, a sum-frequency signal is generated with intensity 
proportional to the cross-correlation of the original pulses as a function of the delay-line 
position. This technique has been used to fully characterize ultrafast pulse propagation*® 
generally disregarding the spatial distribution of the converted signal (cfr. Fig. lb). Nevertheless, 
cross-correlation gating techniques offer several advantages, such as being unaffected by optical 
chirp (in contrast to interferometric methods) and being compatible with collinear schemes 
because possible spurious second harmonic contributions from gate/probe beams, if any, can 
still be spectrally separated. Indeed, a collinear geometry is particularly convenient to mitigate 
possible distortion of upconverted images passing through the BBO, which is a crucial aspect 
for accurate spatially resolved experiments^*. Finally, the gate and probe beams are perfectly 
exchangeable, therefore allowing us to easily investigate any sample at both wavelengths. 
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Upconversion imaging To study the dynamics of light emerging from a complex specimen, 
we image its exit surface on the BBO crystal. There, the original image at the probe wavelength 
is upconverted to a different (signal) wavelength by sum-frequency generation with a gate pulse. 
Non-linear upconversion is a polarization- and wavevector-sensitive process, and both of these 
aspects must be taken into account when designing an imaging apparatus. In the current study, 
which addresses multiple scattering systems, the polarization sensitivity is not relevant because 
all polarization channels have the same optical properties. However, the setup inherently allows 
to investigate separately different polarization channels. In the phase-matching condition, the 
non-linear crystal converts only those wavevectors laying within a relatively narrow paraxial 
range in reciprocal space, which defines the angular acceptance of the BBO crystal (see SI). 
This acceptance angle is weighted differently in the x- and y-directions of the BBO surface 
(that is, along the image plane, see Figs. Ic and d). A double telecentric apparatus in a 4/ 
conhguration ensures that the paraxial component emerging from the sample is preserved at the 
crystal interface. This allows us to obtain a quantitatively accurate reproduction of the original 
intensity distribution at each time step for each point of the object field. This accuracy condition 
holds as long as the transmitted wavevector distribution, however complex, does not depend on 
the point and time of exit, as is the case with multiple scattered light in disordered samples. In 
this case, the intensity carried by any (fixed) subset of wavevectors can be considered to be 
strictly proportional to the total transmitted intensity. 

As seen in Figure Ic, the reduced set of available wavevectors degrades the sharpness of 
the retrieved image differently along the two axes. To obtain the same resolving power in 
both the X and y directions, the aperture of a Fourier filter (FF) can be reduced to act as 
a tunable isotropic wavevector filter. Several alternative upconverting configurations have 
been reported in the literature regarding full-frame imaging applications^’^^^^^, most of which 
employ a non-linear crystal in the Fourier plane of the imaging optics rather than in the image 
plane. Both configurations present different advantages, and the preference for their adoption 
depends on the specific application. In the common Fourier plane configuration, angular 
acceptance limits upconversion uniformity over the entire field of view; in our case, the image 
is formed with a cut-off set of k vectors (thus with lower resolution) while guaranteeing spatially 
uniform upconversion efficiency. Therefore, we adopted the latter configuration because our 
diffuse signal typically does not exhibit any sharp features, and a spatially uniform response 
represents a critical condition for attaining our quantitative accuracy goal over multiple orders 
of magnitude. The generated signal is collected by a lens inside a shielded box where a flip 
mirror allows us to switch between different detection schemes. 

We initially performed standard time-resolved measurements using a photomultiplier tube 
(PMT), that is, collecting intensity without any spatial resolution. Figure lb shows a typical 
cross-correlation measurement of our probe pulse passing through a test target together with a 
fit obtained from the convolution of two squared-hyperbolic secant functions (see SI). . The 
measurement does not exhibit any relevant satellites over several orders of magnitude. By 
switching the photon counter with a CCD camera, a transient image is detected (Fig. Ic) that 
shows the spatial resolution of our setup (~ 22.61pmm“^ in x and ~ 11.31pmm“^ in y at 
0.5 contrast), which is comparable to previous publications^’^^’^^’^^. The final resolution is 
obtained by reducing the aperture of the FF to < 1 mm, until the same (isotropic) resolution of 
~ 11.3 Ip mm“^ is obtained in both directions. 
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With regards to the overall sensitivity of the upeonversion proeess, previous measurements 
performed with the same setup^*’ give an estimate of the signal attenuation at whieh the noise 
level is eventually reaehed. Starting from a typieal pulse energy of ~ 1 nJ at a probe wavelength 
of 810 nm, an upeonverted signal ean still be deteeted after a damping of 8 deeades, taking into 
aeeount the fraetion of the probe beam that is usually lost by diffuse refleetion and the limited 
solid angle subtended by the eolleetion opties. In addition, several other parameters ean be 
adjusted to further enhanee sensitivity, sueh as increasing the integration time of the detector, 
using a non-linear crystal with a higher upeonversion efficiency and a larger angular acceptance 
(for example, bismuth triborate) or increasing the fluence on the crystal. To this purpose, it is 
particularly convenient to increase the gate beam intensity as much as possible to avoid any 
alteration of the sample. 

Data evaluation models By accessing both the spatial and temporal evolution of light emerg¬ 
ing from a complex environment, an incredibly rich picture emerges that needs to be analyzed 
with the proper evaluation tools. Exact modeling of energy transport in disordered samples is 
provided by the radiative transfer equation, which models light propagation as a random walk 
of energy packets undergoing scattering and absorption events. The characteristic length, /g, 
over which a packet is scattered on average is called the scattering mean free path, whereas 
the anisotropic scattering factor, g, expresses the average degree of directionality ranging 
from completely random (g = 0) to perfectly forward (g = 1). Both these parameters are 
directly linked to the microscopic structure of the sample and to the properties of its basic 
scattering elements. However, a similarity relation exists^®’^^ equating a random walk with 
a certain Is and g to a completely isotropic random walk with a rescaled mean free path of 
k = h/i^ - g)- ■ Thanks to this relationship, light transport is often studied in the isotropic 
diffusion approximation (DA) regime, where a single diffusion coefficient, D = kv/S, is used 
as the only transport parameter, with a consequent loss of insight concerning the microscopic 
properties of the sample. Moreover, the convenient description offered by the DA is known 
to be valid only when the transport mean free path is much smaller than the sample size. 
Conversely, its applicability to optically thin media is often debated and a standard approach 
to measure transport parameters in this regime is still missing. Several attempts have been 
made to derive extensions of the diffusive model for optically thin media^^^^^ , and alternative 
techniques to measure transport parameters have been proposed that avoid using diffusion 
theory in the first place^^^^ and are based on single-domain (for example, temporal, spatial or 
angular) investigations. 

In this study, we compare analyses of our experimental data with the mainstream diffusion 
theory and a Monte Carlo model of light transport. While the former method retains its appeal 
due to its simplicity, the latter represents a direct implementation of the radiative transfer 
theory, providing a computationally expensive but asymptotically exact solution to the transport 
problem. This allows us to unveil some of the particularly deceptive ways that the DA can fail 
when dealing with certain types of media. 
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Results and Discussion 

In the following section, we present a few relevant experimental cases where exploiting our 
newly developed ultrafast imaging (UFI) setup, we were able to reveal structural and optical 
features of the investigated samples that could not have been probed unambiguously with 
current state-of-the-art time-resolved techniques and data evaluation models. These claims are 
further validated and quantified by a direct comparison of the experimental data with Monte 
Carlo simulations. The main advantages offered by our experimental technique, as well as 
other possible applications, are discussed in the final subsection. 

Limited thickness artifacts We first tested our setup with the simplest possible disordered 
slab, a homogeneous isotropic sample made of Ti 02 nanoparticles embedded in a polymer 
matrix (see SI). The sample thickness was 203 qm with an average refractive index at SlOnm 
of 1.52. 

Figure 2a shows a normalized time-resolved transmission profile, in addition to a fit using 
the DA relation, which shows excellent agreement with the experimental data. The fit was 
performed with the diffusion coefficient, D, and absorption coefficient, yUa = the only free 

parameters, assuming D = /tv/3 and v = c/n. The fit returns U = 24.1 qm and 4 = 12.0 mm 
(albedo 0.998). The ratio between the thickness of the slab and its transport mean free path 
(referred to as the optical thickness) is commonly used to express the amount of multiple 
scattering in a sample. According to the fit, we obtained an optical thickness of ~ 8.4; a value 
of 8 is usually considered to be the rule-of-thumb lower limit above which diffusive transport is 
regarded as an acceptable approximation"^"^. Nonetheless, the retrieved coefficients predicted a 
total absorption of A ~ 8 %which is much higher than expected because both the host polymer 
and the Ti 02 nanoparticles are known to have vanishing absorption at near infrared wavelengths. 
Even considering a non-absorbing sample with an equal thickness and refractive index, the 
lowest possible decay time that diffusion theory could predict (t° 5 = 6.23 ps, see Fig. 2a, inset) 
would still be appreciably longer than that determined experimentally (t = (6.01 + 0.08) ps). 
The DA accommodates for such discrepancies by introducing an absorption artifact. This can 
intuitively be seen in its relationship for the asymptotic decay time, = n^D/L^^+fx^v (where 
Lgff is the effective thickness depending on the refractive index boundary conditions), the value 
of which can be made arbitrarily small by assuming the presence of absorption. 

In this case, the analysis of the time evolution of the mean square width of the trans¬ 
mitted profile is particularly interesting in that it is inherently free from this absorption-to- 
scattering crosstalk effect^^’"^^ "^^. The mean square width is defined as the variance, w^(t) = 
f r^I{r, t) d^r/ f I(r, t) d^r, where r is a two-dimensional vector in the transmission plane, and 
can be determined if the full spatio-temporal evolution of the profile /(r, t) is known. Due to its 
normalization, any amplitude factor (such as absorption) applied to the profile will cancel out 
exactly at any time, leaving the mean square width unaffected. We therefore further investigated 
the sample using our imaging configuration (which could in principle also substitute the photon 
counter following frame integration, see Fig. 2b). The mean square width, as extracted by the 
obtained frames, clearly shows a linear increase (Figures 2c-d) that can easily be interpreted 
within the DA, which predicts w^if) = ADt. The experimental slope, as obtained from a linear 
fit, corresponds to Df'p = 1746qm^ps“^ k = 26.6 qm, which is appreciably larger than 
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Figure 2: Homogeneous sample data analysis, a Experimental data (black squares) acquired 
with the photon counter (PMT) compared to an exponential fit (solid green line), and a diffusion 
theory fit with D and Ha as free parameters (solid red line). The inset shows the asymptotic decay 
time expected according to the DA for a sample with L = 203 pm, n = 1.52 and //a = 0 
(solid black line). The experimentally retrieved value (dashed green line) is incompatible with the 
DA unless some absorption is introduced, b A Monte Carlo (MC) solution (solid red line) of the 
radiative transfer equation showing perfect agreement with the experimental data for a non-absorbing 
sample. The integrated intensity from UFI frames and the PMT data (cyan circles) are shown, c A 
set of frames acquired at different delays with our UFI setup. Each frame is averaged over different 
disorder realizations and is displayed normalized to its own maximum intensity, d Comparison of the 
experimental and simulated mean square width. 
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the value retrieved from the time-resolved eurve. We therefore ran Monte Carlo simulations 
(see SI) to resolve this diserepaney, and eventually sueeeeded in simultaneously reprodueing 
both the experimental time-resolved and mean square width evolution datasets using a single 
transport parameter, 4 (Figures 2b and d). The Monte Carlo inversion proeedure gives a value of 
k = 25.5 qm (again eorresponding to an optieal thiekness of ~ 8) whieh is perfeetly eompatible 
with both datasets without introdueing absorption. 

It should be noted that the transport mean free path, as retrieved from the mean square width 
slope in the diffusive framework, results in just a 4 % overestimation of the aetual best value 
retrieved by the Monte Carlo simulations. Sueh a level of aeeuraey is aeeeptable for a wide 
number of applieations, where UFI teehniques eould be exploited to extend the applieability of 
the DA in this intermediate thiekness range. Moreover, as widely suggested in the literature"^^’"^^, 
estimating transport properties from the mean square width rather than from time-resolved data 
is mueh more aeeurate and straightforward beeause in the former ease, a priori knowledge of 
the absorption eoeffieient, refraetive index eontrast and sample thiekness is not required. Indeed, 
the simple ADt dependenee does not involve these parameters or the extrapolated boundary 
eonditions. This simple determination of the slope is also insensitive to the exaet determination 
of the time axis origin and the exaet time-window eonsidered for the fitting, both of whieh 
represent long-standing issues in the evaluation of time-resolved data"^’’^^. In addition, for the 
same reason that the mean square width is independent of the absorption, it is also effeetively 
independent of the integration time, as well as of any laser power fluetuations. Interestingly, 
the time-resolved DA ean deeeptively eompensate for its limited aeeuraey by introdueing an 
absorption artifaet, but laeks a solution when the eorreet value of absorption is eonsidered. This 
point is partieularly relevant for investigations of optieally thin speeimens, sueh as biologieal 
tissues"^^. 

As previously mentioned, our UFI teehnique is eapable of replaeing eompletely integrated 
measurements with eomparable sensitivity (Fig. 2b). Exploiting the full spatio-temporal 
information, the resulting set of frames ean be deeomposed into a mean square width and a 
time-resolved eurve, providing two almost independent datasets in a single measurement session. 
Indeed, the former does not depend on the integrated intensity of the individual frames, whereas 
the latter is primarily linked to the transport along the z axis. In this respeet, our perfeetly 
matehing Monte Carlo validation, whieh is eondueted independently on two experimental 
datasets that depend differently on different sets of parameters, is a elear indieation of the 
quantitative aeeuraey and reliability of our setup over many orders of magnitude in both the 
temporal and spatial domains. 

Large-scale heterogeneity artifacts We also used our UFI teehnique to investigate more 
diverse and eomplex media. A new polymer slab (L = 190 qm) with a roughly doubled 
seatterer density was manufaetured to obtain an optieally thieker sample. Again, a perfeet fit 
eould be obtained in the diffusive framework (Fig. 3a) if the presenee of absorption was ineluded 
at an unlikely integrated value for the total absorption, A ~ 7 % (/t = 12 qm, 4 = 13.2 mm). 
In sharp eontrast to the previous experiment, a non-absorbing Monte Carlo fit was unable 
to reproduee the time-resolved data, even when the observed experimental deeay time was 
perfeetly matehed (Fig. 3b, dotted blaek line). Furthermore, the mean square width exhibited a 
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Figure 3: Heterogeneous sample data analysis, a Photon counter (PMT) data (black squares) 
and best diffusive (DA) fit (solid red line) giving /( = 12 gm and 4 = 13.2 mm. b, c Comparison 
of the PMT and UFI data to a multi-layered simulation (solid red line) and two homogeneous 
simulations. The black-dotted line refers to the simulation that reproduced the experimental decay 
time (/, = 14.2 gm), while the black-dashed line shows the simulation that reproduced the mean 
square width slope (/, = 15.7 gm). The decay-time-reproducing simulation fails to reproduce the 
mean square width slope, while the simulation that shows perfect agreement with the mean square 
width slope poorly matches the time-resolved decay, therefore revealing the inconsistency of the 
homogeneity hypothesis. The Monte Carlo (MC) curves were rescaled to match the experimental 
peak count values. 


steeper than expeeted slope (£>exp = 1050 gm^ ps“^ ^ /j = 16 gm) as if the diffusion proeess 
was enhaneed along the in-plane direetions, as shown in Figure 3e. 

Driven by this diserepaney, we further investigated the sample with a seanning eleetron 
mieroseope (SEM), revealing a heavily layered modulation of the seatterer density eompared 
to the more homogeneously dispersed one in our first sample (Figures 4a and b). Taking 
advantage of the versatility of Monte Carlo simulations, we modeled our sample after our SEM 
images as being eomposed of hve layers with three different thieknesses and densities arranged 
symmetrieally with respeet to the eentral layer. A new Monte Carlo brute foree ht assuming a 
hxed geometry and three free parameters, ki, la, k-i, was eventually able to perfeetly reproduee 
both the time-resolved and mean square width data simultaneously, as shown in Figures 3b- 
e. We obtained transport mean free paths of In = 3.5 gm, /,2 = 21.5 gm and = 11 gm, 
respeetively, for the high, low and intermediate-density layers. The absorption eoeffieients were 
all set to //a = 0 gm“4 In this ease, the mean square width/time-resolved deeomposition elearly 
allowed an unexpeeted degree of eomplexity to be revealed and eharaeterized thanks to the 
rigid eonstraints of the double boundary datasets. Resolving this eomposite strueture showed 
that a layered heterogeneity ean, in prineiple, mimie the effeets of anisotropie transport and 
therefore revealed a previously unaddressed physieal effeet. This is espeeially relevant given 
the pervasiveness of layered media (for example, eoatings, atmospherie physies and biologieal 
tissues), whieh often exhibit eounter-intuitive features that standard time-resolved teehniques 
are still unable to explain, as in the ease of light transport aeross the human forehead^'. 

Ultrathin samples Despite the advantages offered by deeomposition into a mean square width 
and a time-resolved eurve, it is elear that the raw, non-integrated data provide an irredueible set 
of information. A dramatie illustration of this point is obtained by studying a thin biologieal 
sample, sueh as a small strip eut from the dried skin of a sliee of grapefruit. As the SEM image 


9 



























b) homogeneously disordered sample 



c) grapefruit slice dried skin 


Figure 4: SEM side-view images of the anaiyzed sampies. a Heterogeneously disordered sam¬ 
ple. The central inset shows a transverse section, which exhibits a layered symmetric structure 
that we modeled with a high-density central layer (Li = 26 pm), a low density interstitial layer 
(L 2 = 60 pm) and an intermediate density outer layer (L 3 = 22 pm). The heterogeneity represented 
by the core layer is also appreciable at the optical microscope (left inset). The right insets highlight 
the density difference between regions 1 and 2. b A transverse section of the homogeneously 
disordered reference sample. Despite showing signatures of scatterer clustering, small-scale clus¬ 
ters/heterogeneities are homogeneously dispersed across the sample thickness and do not exhibit 
any layered arrangement. As confirmed by the excellent agreement obtained with a single parameter 
Monte Carlo simulation and other work on similar samples®”, light perceives an effective, homoge¬ 
neous scatterer density in contrast to the previous tightly layered sample, c Lateral section of the 
dried skin of the grapefruit slice. 


reveals (Figure 4e), its structure consists of a conglomerate of small flakes forming a corrugated 
slab ~ 85 pm thick. The dried sample appears to be a brittle, almost transparent membrane 
(Figure 5a, inset). 

In this range of optical (and absolute) thicknesses, standard experimental techniques fall 
short because of the extreme time scales involved, causing common diffusive modeling to 
fail dramatically. Furthermore, most of the signal will be ballistically transmitted, carrying 
almost no information concerning the optical properties of the sample. As Figure 5a shows, the 
light that is scattered inside the grapefruit membrane is l/50th of the original intensity. Our 
ultrafast imaging technique allows us to selectively address the light that was retained for a 
longer time inside the sample whereas at the same time spatially inspecting it as it propagates 
along the main (transverse) slab extension. Figure 5b shows a collection of frames acquired 
over a time-window of ~ 5 ps after pulse injection, corresponding to a total path length greater 
than 13 times the sample thickness. We see that the light spreads through the membrane with a 
well-dehned wavefront traveling inside its main plane, resulting in a dramatic departure from 
any prediction compatible with the diffusive framework. However, standard time-resolved or 
steady state investigations would still just measure a single decay time and a bell-shaped prohle, 
respectively, which could deceptively support an inappropriate interpretation in terms of the 
DA. 

As we will show, our ultrafast time- and space-resolved experiment provides a deep insight 
into the light transport properties of complex systems, which requires the development of a 
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Figure 5: Grapefruit skin data anaiysis. a The time shift induced by the sample in the probe path 
(black circles vs. black squares). Knowing its thickness from previous scanning electron microscopy 
images, the effective refractive index is found as r = 1 + cAtjL ~ 1.31. The inset shows the actual 
appearance of the investigated specimens, b The frame set, each of which is averaged over a small 
portion of the sample. The transmitted intensity distribution shows a quasi-ballistic elliptic wavefront 
traveling at slightly different speeds (see SI), c Horizontal cross-cuts of the wavefront and a Monte 
Carlo fit with g and 4 as free parameters. The fit returns values of ^ = 0.7 and 4 = 150 pm. d 
Simulated radial intensity profiles in this extremely optically thin regime show a significant breakdown 
of the similarity relation. 
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novel analysis methodology. We demonstrate a proof-of-eoneept proeedure to assess in-plane 
transport properties in this extremely optically thin regime. For the sake of simplicity, we 
will illustrate our investigation on a cross-cut of the profile assuming isotropic transport. The 
general treatment is analogous and will be the subject of future work. Figure 5c shows a 
comparison between experimental cross-cuts at different times and the corresponding time 
evolution of an isotropic Monte Carlo simulation that succeeds in reproducing the wavefront 
speed, intensity decay and overall shape at all given times. To achieve this collective agreement, 
the anisotropic scattering coefficient, g, needs to be considered as an additional parameter 
to express the degree of directional correlation between consecutive scattering events. As 
mentioned, diffusion theory is degenerate in g as long as the transport mean free path 4 is 
kept constant, as expressed by the similarity relation k = 4/(1 - g)- This is why retrieving 
any information concerning the angular scattering function is particularly challenging, and 
one must therefore resort to a regime where the DA ceases to hold"^^. Indeed, as shown by 
Figure 5d, the observed transverse intensity patterns exhibit a significant breakdown in this 
degeneracy condition in thin slabs, with g playing a substantial role in determining the shape 
and evolution of the traveling wavefront, even when 4 is constant. In particular, while the outer 
wavefront propagates almost ballistically (see SI), its instantaneous position and peak intensity 
vary appreciably with different combinations of g and 1^, therefore allowing both to be retrieved 
with an estimated precision of ~ 10%. The set of simulated curves shown in Figure 5c is 
obtained for g = 0.7 and 4 = 150 p,m (corresponding to an in-plane 4 = 500 p,m). However, 
quantitatively accurate figures for this particular specimen would require an analysis involving 
fully anisotropic simulations. 

The main purpose of this experiment is to showcase the vast possibilities opened by rich 
time- and space-resolved outputs. Our g estimation results from a collective fitting routine 
involving multiple spatio-temporal datasets and stands in contrast to common and less robust 
g determination techniques involving a single scalar, such as the attenuation coefficient of 
a collimated beam^^’^^. Accurate determination of single scattering anisotropy is crucial to 
directly probe the shape, orientation and optical properties of single microscopic scatterers and 
might therefore be of interest to a broad range of fields, from realistic computer rendering of 
participating media^"^’^^ to photo-therapeutic diagnostics of biological tissues^^’^^. Moreover, 
compared to other available techniques, the determination of g and /g is achieved in a single 
measurement session, which also simplifies the retrieval procedure. The peculiar doughnut¬ 
shaped intensity pattern that enables such a rich characterization of the scattering properties 
is of course not unique to this particular sample, and we observed it in a variety of different 
biological samples of both plant and animal origin, usually differing only by the actual contour 
shape of the propagating wavefront. This supports the general validity and usefulness of 
our proof-of-concept evaluation technique in studying extremely complex media, included 
biological tissues. 

Discussion As illustrated, our optically gated imaging technique offers several advantages, 
including higher quantitative fidelity over a large field of view with respect to Fourier plane 
configurations and sub-ps time resolution compared to other electronically gated imaging solu¬ 
tions. Gaining direct access to such time scales not only enables the study of optical properties 
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in thin or inherently minute specimens (for example, ocular fundus, vascular walls, skin dermis 
and dental enamel), but is fundamental to investigate light transport at the mesoscopic level. 
From a technical viewpoint, the main strength of our setup is its simplicity and lack of complex 
electronics. The inherently synchronous nature of the pulses renders the entire setup insensitive 
to timing jitter. Power fluctuation/drifts also do not pose a problem when measuring the mean 
square width evolution, which also allows us to dynamically adjust the integration time to 
enhance the signal-to-noise ratio. This flexibility should make it possible to replicate the setup 
with cheaper and more easy to use hands-off fiber lasers, therefore widening the effective 
applicability of the technique. 

Accessing both transverse and axial transport simultaneously through a full spatio-temporal 
investigation allows us to correctly evaluate the applicability of the common diffusion frame¬ 
work, which can be subject to appreciable absorption overestimation in samples presenting a 
limited optical thickness or large-scale heterogeneities. An accurate retrieval of the absorption 
coefficient is particularly relevant in many biomedical applications where its value is directly 
linked to quantities such as chromophores concentration or blood oxygenation levels. The 
observed discrepancy, which would be not apparent with state-of-the-art spatial or temporal 
measurements, is revealed by comparing standard time-resolved transmitted intensity decay 
to the evolution of the mean square width, both of which can be acquired simultaneously. 
Nevertheless, the most relevant advantage offered by direct access to the transient mean square 
width expansion is that it enables a robust and simple interpretation. In sharp contrast to other 
common observables, such as the total transmittance and its decay time, the mean square width 
is independent of absorption and largely unaffected by refractive index contrasts and sample 
thickness, therefore overcoming long-standing problems posed by their precise assessment. 
One prominent case where this feature could be decisive is that of the experimentally observed 
decrease in the diffusive coefficient of a turbid slab with decreasing thickness^*, which was 
indirectly determined through the transmittance decay time. A spatio-temporal measurement 
in terms of the mean square width evolution would allow us to directly probe the diffusion 
coefficient irrespective of the exact boundary conditions, the incorrect estimation of which has 
been cited as a possible cause for this apparently anomalous behavior^^. 

A different scenario emerges when studying samples with extremely low optical thickness, 
which revealed a peculiar transport regime where /g and g, rather than 4, are the meaningful 
parameters determining how radiative energy spreads through the sample. This finding illus¬ 
trates a further practical situation where light transport must be studied as a spatio-temporal 
phenomenon in its entirety. This point has not yet been properly considered when studying 
transport in thin samples. Moreover, we have demonstrated how this transport regime enables 
accurate retrieval of otherwise barely accessible transport parameters, especially for this class 
of samples where well-established experimental techniques are still missing. 

In general, several applications will be facilitated by the three-dimensional characteriza¬ 
tion capabilities of the setup (see also SI) opening new possibilities to properly investigate 
structurally anisotropic media that are of paramount relevance in many industrial and biomed¬ 
ical applications, and whose correct theoretical modeling in terms of light transport is still 
an object of lively debate in the literature^^’^®. Other relevant open questions that might be 
tackled directly with a spatio-temporal investigation include experimental verification of the 
polarization dependence predicted for the diffusion coefficient^^ as well as the occurrence of 
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non-monotonic mean square width evolution in samples exhibiting Anderson localization"^^ 
or the lateral expansion of Levy type transport, where the experimental limit posed by the 
finite thickness of the sample might be circumvented^^. Similarly, an UFI measurement of 
intensity profiles reflected by macro-porous bulk materials could access light spreading within 
the first few hundred femtoseconds from which the time-dependent diffusion constant could 
be determined^^ leading to information concerning the structure factor of the material, which 
would not be obtainable at long times. 

Conclusions 

The field of light transport in complex or structured media still presents a multitude of ex¬ 
perimental challenges that could significantly benefit from the advent of novel multi-domain 
investigation techniques. In this study, we demonstrated an experimental technique to go be¬ 
yond the ps resolution barrier for spatial investigations of transient intensity profiles emerging 
from arbitrary samples. Gaining access to this time scale with a broad field of view is key to 
the characterization of a large class of samples, ranging from thin biological membranes to 
complex optical microdevices. As we demonstrated, spatio-temporal investigations enable the 
study of heterogeneous or semi-transparent media that are difficult to characterize with present 
state-of-the-art techniques. Further, when dealing with more standard samples, combined 
access to both transverse and axial transport helps reveal the presence of flaws or artifacts in the 
common diffusion-based data evaluation framework. Remarkably, in many practical situations, 
we found that low optical thickness can be exploited as an advantage rather than a drawback, 
therefore offering a complementary framework to study the broad class of samples that do not 
meet the validity conditions required by the DA. Nonetheless, the setup that we propose is 
generally applicable and non-invasive, usable at different wavelengths, with a broad field of 
view and high-time resolution. We therefore envision that it could be applied to a wide range 
of photonic applications for both the sub-ps physics it allows to investigate and its convenient 
wide-field acquisition, which does not require scanning over the region of interest. 
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Supplementary Information to: Spatio-temporal visualization of 
light transport in complex photonic structures 

Experimental setup details 

The non-linear erystal used is a square 5 mm x 5 mm x 2 mm BBO erystal. Foeal lengths 
are respeetively 100 mm for L2 and L3, and 125 mm for the lens eolleeting the up-eonverted 
signal after the BBO. The final magnifieation obtained through the double imaging stages was 
roughly lx. The foeusing lens LI is ehosen so to provide a point-like exeitation spot with 
respeet to the typieal length seales of the sample, whieh in our ease was set to ~ 10 qm. As 
far as the mentioned eonstraints on the waveveetor distribution manipulation are satisfied, a 
different set of opties ean be employed to obtain a different magnifieation. The CCD eamera is 
a baek-illuminated Andor iKon M912, with a 512x512 pixel sensor. 

Polymer samples fabrication 

Investigated polymer samples are made of a commercial UV-curing acrylate optical adhesive 
(Norland 65) with a dispersion of rutile nanoparticles with a diameter of 280 nm (Huntsman’s 
Tioxide R-XL). The mixture of polymer and nanoparticles is rendered homogeneous through 
magnetic stirring and a ultrasonic bath (~ 1 h). We then let the resulting opaque paste be sucked 
by capillary forces inside an air gap of controlled thickness formed between two microscope 
glasses. The two glasses are firmly held apart by micro-spherical spacers of calibrated size. By 
spin coating in advance the microscope glasses with water-soluble polyvinyl alcohol we are 
able at a later stage to remove the glasses obtaining a flexible, free standing polymer slab. This 
notably allows to avoid multiple reflection of both pump and scattered light coming from the 
enclosing glasses. 

Monte Carlo software 

All simulations shown have been performed with a C-I--I- software called MCPlusPlus^ based 
on the standard MCML Monte Carlo code for multi-layered slab samples^. Its source code 
is currently under development and is freely available at http://www.lens.unifi.it/ 
quantum-nanophotonics/mcplusplus/. Benefiting from the object-oriented paradigm 
of C-I--I-, the software offers a high level of abstraction, scalability, modularity and ease of 
maintenance. Being completely written in C-I--I-, it can be executed an any hardware and can 
take extensive advantage of the multi-thread capabilities of most modern CPUs for increased 
performance. Other features worth mentioning include python scriptability and both raw and 
histogrammed output in the convenient H5 binary format. 

Probe and gate pulses characterization 

A Spectra Physics Tsunami Ti:Sa laser is used to pump a Opal optical parametric oscillator, 
which provides a signal at (1510+ 11) nm, an idler (not used in the experiment) and an 
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Figure S1 : Pulses characterization, a Spectrum of the unconverted residual (probe pulse) from 
the OPO, with a sech^(a») fit. b Spectrum of the OPO signal (gate pulse) with a sech^(aj) fit. c 
Autocorrelation measurement of the residual pulse. The envelope profile fringes has been fitted with 
the autoconvolution of two sech^(f) pulses. 


unconverted residual at the same wavelength of the pump (i.e. (810 + 5)nm). The average 
power of the signal and the residual beam is respectively 190mW and 280 mW. Both pulses 
are characterized at the beginning of each measurement session spectrally and temporally. 
Figures Sla and b show typical spectra recorded for the Tsunami and the Opal, fitted with a 
I{oj) oc sech^(|(a» - 6 Uc)Tp) function, as expected for an actively mode-locked Ti:Sa laser in a 
sub-ps configuration (i.e. featuring ODD compensation). Figure Sic shows an autocorrelation 
measurement of the residual pulse, fitted with the autoconvolution of two sech^ pulses, which 
can be expressed as /ac (0 = 7 COsech^( 7 t) [ 7 COtanh( 7 r) - 1 ] with 7 = ( 21 og(l -l- V 2 ))/(Tfwhm)- 
The full width half maximum duration of the original pulse is obtained by multiplication for the 
appropriate deconvolution factor, giving Tfwhm = 0.6482 x 144.6 fs = 93.7 fs. To estimate the 
duration of the gate pulse we can rely on the cross-correlation measurement. Due to the lack of 
any analytical expression for the cross-convolution of two different sech^ pulses, we performed 
a fit plugging at each iteration the numerical convolution of the known pump pulse with an 
unknown sech^ pulse (Fig. lb). The routine eventually returned a full width duration of 134 fs 
for the gate pulse. The cross-correlation was measured with and without the lenses, to make 
sure that the optics used for the experiment were introducing negligible dispersion. The final 
full width half maximum, representing the instrument response function is of approximately 
170 fs. The source term in all Monte Carlo simulation showed in this paper has been modeled 
as a sech pulse of this duration, while its spatial distribution was set to a Gaussian beam with 
a waist of 10 ^im, as given by our focusing lens. 

Full 3D retrieval of ballistic speed of light 

Spatio-temporal measurements of light traveling inside the grapefruit membrane and other 
analogue samples exhibit a transmitted intensity pattern which is determined by both ballistic 
and scattered light — where by “ballistic” we refer to light that undergoes roughly just two 
main scattering events: one to get inside the membrane and one to be scattered out. In particular, 
such ballistic component appears to persist for the whole time scale observed in our experiment 
and can be easily addressed by measuring the instantaneous position of the outer wavefront of 
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Figure S2: Grapefruit in-piane iight speed retrievai. a For each measured profile, the instanta¬ 
neous position of its outer wavefront is calculated at half Its maximum amplitude A for both horizontal 
(shown) and vertical cross-cuts, b For the whole time window taken into account by our mea¬ 
surements the outer wavefront appears to travel ballistically with a well-defined speed. Linear fits 
return speeds of 170 pm ps’^ (0.57c) and 200 pm ps“^ (0.67c) for the horizontal and vertical axis 
respectively. 


the traveling pulse where, by definition, it is aeeumulated. This is eonfirmed by the steadily 
linear inerease of the wavefront position. The retrieved slopes direetly give the light speed along 
a eertain direetion inside the slab speeimen, whieh notably appears to be different along the x 
and y direetion for the investigated sample. Combining these results with the perpendieular 
delay introdueed by the intervening sample whose thiekness is known (efr. Figure 5a) we are 
able to retrieve the speed of light inside the membrane along eaeh spatial dimension. 
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